Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 6 August 2012 (MN KTeX style file v2.2) 



Radial Migration in Disk Galaxies I: Transient Spiral Structure and 
Dynamics 



(N 

o 

bJQi 



o 

6 

5— ( 

(N 
> 
m 



Rok Roskar 12 *, Victor R Debattista 3 ' 4 , Thomas R. Quinn 2 , James Wadsley 5 

1 Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland 
2 Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA 
3 Jeremiah Horrocks Institute, University of Central Lancashire, Preston, PR1 2HE, UK 
A RCUK Fellow 

5 Department of Physics and Astronomy, McMaster University, Hamilton, ON, L8S 4M1, Canada 



6 August 2012 



ABSTRACT 

We seek to understand the origin of radial migration in spiral galaxies by analyzing in 
detail the structure and evolution of an idealized, isolated galactic disk. To understand the 
redistribution of stars, we characterize the time-evolution of properties of spirals that sponta- 
neously form in the disk. Our models unambiguously show that in such disks, single spirals 
are unlikely, but that a number of transient patterns may coexist in the disk. However, we also 
show that while spirals are transient in amplitude, at any given time the disk favors patterns 
of certain pattern speeds. Using several runs with different numerical parameters we show 
that the properties of spirals that occur spontaneously in the disk do not sensitively depend 
on resolution. The existence of multiple transient patterns has large implications for the or- 
bits of stars in the disk, and we therefore examine the resonant scattering mechanisms that 
profoundly alter angular mo menta of individual stars. W e confirm that the corotation scatter- 
ing mechanism described by Sellwood & Bi nnevl d2002l) is responsible for the largest angular 
momentum changes in our simulations. 
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1 INTRODUCTION 

As stellar galactic disks form and evolve, the processes govern- 
ing their development leave behind a multitude of traces. Massive 
stars pollute the interstellar medium through supernova explosions 
and stellar winds, endowing subsequent generations of stars with 
distinct chemical signatures. Similarly, dynamical perturbations, 
whether they are secular (spiral arms or bars) or due to the hier- 
chical build-up of mass, sculpt the kinematic properties of stellar 
populations. Together, these signatures provide a "fossil record" of 
a disk's evolution through cosmic time. 

Galaxy formation models attempt to match their predictions to 
this fossil record and attempt to reconstruct the disk's history. For 
the past thirty years, an implicit assumption has been made in such 
modeling: star s remain in the same part of t he disk forever (e.g. 
Tinslev| |l975|: |Matte ucci & Francoisl 1 19891 : iBoissier & Pr antzos 
19991 : IChiappini et alJl20oC This simple assumption carries enor- 
mous power as it allows one to reconstruct the time evolution of 
a given quantity (such as metallicity and star formation rate) at 
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a particular radius from present-day, single epoch observations. 
Such "static" modeling has been successful in illuminating sev- 
eral critical aspects of disk formation such as the n eed for infall o f 
low metallicity gas to solve the G-dwarf problem ( iTinslevll 197^1 . 
Over the past several years, however, this assumption of stars re- 
maining near their birth radii has been firmly shaken by the re- 
alization that rapid stellar migrations of several kpc are possible 
dSellwood & Binnevll2002l. SB02 hereafter). 

The idea that stars may not remain near their birth radii is not 
new. IWielenl d 19771) suggested that stellar orbits diffuse in veloc- 
ity space consequently inducing a drift in the gal actocentric radius 
for an ensemble of stars. This was followed up bv lWielen & Fuchisl 
(1985) and applied to the determinat ion of the Sun's birthplace 
dWielen et alJI 19961) . In lWielenl dl977l) the diffusion is a relatively 
slow process driven by random scattering by giant molecular clouds 
- for observationally-constrained velocity dispersions, variations in 
galactocentric radius of at most few kpc are expected. Instead, in 
the mechanism described by SB02, migrations of several kpc may 
take place in a few hundred Myr due to very efficient exchange of 
angular momentum at the corotation resonance (CR). 

This exchange of energy and angular momentum at the CR 
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occurs without changing the orbital radial actions therefore retain- 
ing the orbital circularity. In addition, the mechanism operates most 
efficiently for stars on the most circular orbits. These two proper- 
ties lead to the peculiar result that migrations of several kpc in a 
MW-size disk are possible without substantial accompanying in- 
crease in random mot ion. The increase in v elocity dispersion with 
age of stars (see e.g. iHolmberg et alj|2009l for the solar neighbor- 
hood) constrains heating processes and thereby also the magnitude 
of spiral perturbations, so the ability of the CR to redistribute stars 
substantially without heating excessively is a critical aspect of the 
CR mixing process. Because the SB02 mechanism involves the CR, 
the spiral amplitudes must be transient if redistribution rather than 
trapping is to occur. In the absence of other perturbations, a steady 
spiral traps stars on horseshoe orbits - however, if a spiral's ampli- 
tude grows and decays on a timescale comparable to one half the 
libration period of a horseshoe orbit, the spiral will merely deposit 
the star on the other side of the resonance and vanish before pulling 
it back (SB02). 

Unlike steady spirals that can heat the disk o nly at 
the Lindblad resonances dLvnden- Bell & Kalnaisl Il972h . tran- 
sient spirals c an also heat the disk away from the princi- 
pal resonances l Barbanis & Woltidll967l : iBinnev & Lacevfl 19881 ; 
Ijenkins & Binnevll990h . Consequently, a single spiral can heat the 
disk enough to prevent any furth er asymmetric str ucture formation 
in just a few rotations (Sell wood & Carlberg 1984). Star formation 
provides a natural cooling mechanism, continually peppering the 
disk with young stars bom with the galactic velocities of their par- 
ent gas clouds. ISellwood & Carlb erg ( 1984) found that in the pres- 
ence of such cooling, transient spirals keep the disk in a quasi-stable 
state at a Toomre Q ~ 2, which allows for their continual regener- 
ation. 

The realization that stars in galactic disks may migrate ra- 
dially across significant distances has in recent years completely 
changed the discourse on spiral galaxy evolution. Following SB02, 
several other theoretical works have further illuminated the com- 
ple xities that rad i al mi xing introduces for stellar population stud- 
ies. iLeoine et alj d2003l) investigated t he effect of corotation scat- 
tering on the disk metallicity gradient. Ros kar et al.1 d2008bl here- 
after R08a) showed that radial migration could drastically alter the 
stellar population properties of outer disks and provided an expla- 
nation for the observed g radients across the s urface brightness pro- 
file break in NGC4244 dde Jong et ai1l2007l) . R08a predicted that 
disks with broken exponential profiles should show an inflection 
in the mean age profile corresponding to the break radius. This 
has subsequently been confirmed indirec tly by surface photometry 
(Bakos et al. 2008]; |AzzoUmirt^]l2008b and directly by integral- 
field spectroscopy dYoachim et alJl2010L Yoachim et al. 2012, in 
press) and resolved-star counts (Radburn-Smith 2012, submitted), 
providing further evidence that radial mixing occurs in external 
galaxies. 

Using the same simulations, Roskar et al. (2008a, hereafter 
R08b) extended the analysis of R08a and investigated the reper- 
cussions of stellar migration for a range of stellar population stud- 
ies: the solar neighborhood age-metallicity relation (AMR) and 
metallicity distribution function (MDF), the evolution of metallic- 
ity gradients, and the reconstruction of star formation histories from 
present-day observations of stellar populations in external galax- 
ies. They found that > 50% of stars on mostly circular orbits in 
the solar neighborhood of their model have come from elsewhere 
and that migration flattens the AMR and broadens the MDF. They 
also found that the reconstruction of a star formation history be- 
comes problematic in the presence of migration especially at large 



radii where the mass in migrated stars approaches or exceeds the 
cumulative amount of stars formed in-situ. R08ab therefore estab- 
lished the notion that regardless of the method used to observe a 
galactic disk (unresolved photometry, photometry of resolved stars, 
spectroscopy of individual stars in the MW) radially migrated stars 
substantially alter the combined properties of the observed sample, 
thereby possibly strongly biasing the outcome of any subsequent 
modelling. 

ISchonrich & Binnevl d2009tJfbh included a probabilistic pre- 
scription of radial migration in a chemical evolution model con- 
strained by Milky Way observables and found that, similarly to 
R08b, the MDF is broadened and the AMR flattened as a result 
of migration. They found also that most obervational properties 
and pec uliarities of the thick disk can be explained by radial mi- 
gration. lLoebman et alj d201 II) explored the formation of the thick 
disk via radial migration in the iV-body models of R08ab and sim- 
ilarly found that simulated trends tend to agree qualitatively with 
observed thick disk propert i es from the Sloan Digital Sky Survey 
(SPSS ) dlvezic et alj 120081; |Lee et alj 1201 ll) . iMinchev & Famaevl 
(2010); Min chev et alj (2011, 2012) argued that apart from corota- 
tion scattering as presented by SB02, combined effects of multiple 
nearby resonances from several different patterns may also be im- 
portant in driving the redistribution of stars. Their results implied 
that subst antial mixing is possi ble even if the spiral structure is not 
transient. iBrunetti et al.1 J20TH) measured the diffusion coefficients 
in idealized bar-unstable disks and found the CR of the bar to be 
driving the diffusion, though their models were collisionless and as 
such did not allow for recurrent spiral activity. 

Understanding radial migration is particularly relevant at the 
present time because of the upcoming surveys designed to study the 
detailed struct ure of the Milky Way disk. S tudies using data from 
the SDSS (e.g. lJuric et ai]|2008l;[lvezic et al|2008l) a nd its follow- 
up surveys such as SEGUE ( Yannv et all 120091 ; Ide Jong et al.l 
l20ld:lLee et alj|201lh. as well as results using data from RAVE 
dWilson et alj 1201 ll ; iRuchti et al.1 1201 ll) have already pushed the 
current models to their limits in trying to explain the various ob- 
served trends and interdependencies in the thick and thin disks. 
In the coming years, spectroscopic surveys, such a s HER MES 
( Free man & Bland-Hawthornl 12008b and APOGEE dPrieto et al.l 
2008), will obtain high-resolution spectroscopy of millions of stars, 
allo wing finally for "chemic al tagging" of stars into their birth clus- 
ters dFreeman & Bland-Hawthorr]|2002l) . Such observations will in 
principle allow for a d irect measurement of stellar radial migra- 
tion in the Milky Way dBland-Hawthorn et alJlioToh . At the same 
time, the Gaia mission will p rovide a complete 6D map of a 10 kpc 
sphere centered on the Sun dPerrvman et alj|200ll ). and follow-up 
spectroscopic surveys will yield vast amounts of complementary 
chemi cal abundance data. Finally, the Large Synoptic Survey Tele- 
scope dLSST Science Collaborationsl l2009) will similarly provide 
ground-based photometry of the entire sky. To distill these vast data 
sets and apply them to galactic archeology, a clear understanding of 
the radial mixing processes is essential. 

In this paper, we explore in detail the properties of self- 
propagating spiral structure and the cause s of the resulting radia l 
migration in simulations from R08ab and lLoebman et alj ( 201 ll) . 
We first quantify the spiral structure that forms spontaneously in 
our simulations and use this analysis as a basis for understanding 
the causes of radial migration. We also perform a set of numerical 
tests to explore the robustness of our results to choices of numerical 
parameters and stochasticity. The Paper is organized as follows: in 
§|2]we discuss the details of our models; in §|3]we quantify the spi- 
ral structure in our fiducial simulation; in §[4] we explore the causes 
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of radial mixing; in §[5] we present several tests of the effects of 
numerical parameters on the generation of spiral structure; finally, 
we state our conclusions in §|6] 



2 SIMULATIONS 

The initial conditions for all of the runs presented here are gen- 
erated as described in iKaufmann e t al. (2006, 120071) . and consist 
of spheri cal distributions of d ark matter (DM) and gas following 
an NFW dNavarro et al.ll 19971) density profile. The random veloc- 
ities of the DM particles are initialized to ensure equilibrium by 
means of the distribution function obt ained from an inversion of 
the NFW density profile dKazantzidis et alj|2004l) . The gas is ini- 
tialized to the same mass distribution with a temperature profile to 
yield an approximate hydrostatic equilibrium. The gas is also given 
a spin consistent wit h values obtained from collisionless c osmo- 
logical simula tions (Bu llo cket al.ll200ll; iMaccio et al]|2007n . i.e. 
A = (J/G) V / \E\/M 5 = 0.039, where J is the total angular mo- 
mentum, E is the total energy of the system and M is the mass. 
The angular momentum follows a radial profile of j tx r, where j 
is the specific angular momentum. In the fiducial run we use 10 6 
particles per component, resulting in DM particle mass of 10 6 M© 
and initial gas mass of 1.4 x 10 Mq. 

The simulation s were run using the code GASOLINE 
dWadslevetai] r2004), a hydrodynamics e xtension of t he parallel 
multi-stepping N-body code PKDGRAV dStadelll200ll) . The tun- 
ing of the sub-grid star formation and feedb ack prescriptions are 
described extensively in lStinson et ~al] d2006l) . Here we summarize 
the important features. A gas particle becomes eligible for star for- 
mation when its density exceeds 0.1 cm -1 and its temperature dips 
below 15,000 K. If eligible, the gas particle converts some of its 
mass into a star particle at a rate given by 

«£ =c ftfi., (1) 

where c* is the star formation efficiency (set to 0.05 in all of the 
simulations), p gas is the gas density and td yn is the local dynami- 
cal time. Star particles form with 1/3 of the initial gas particle mass 
(in the fiducial simulation this translates to star particle masses of 
4.7 x 10 4 Mq) and to avoid unreasonably small particle masses, 
gas particle masses are limited to 1/5 of their initial mass. When a 
particle crosses this threshold its mass is distributed among neigh- 
boring particles, resulting in an overall decrease of the number of 
gas particles with time. Each star particle represents an entire stellar 
population and therefore a spectrum of stellar masses d escribed by 
the Miller-Scalo initial mass function ^Miller & Scalol l 19791) . The 
evolution of massive stars is followed and a feedback cycle is ini- 
tiated to reflect the explosion of type II supernovae. At typical par- 
ticle masses of ~ 1O 4 M0, the supernova energy is injected into 
the ISM "in bulk", i.e. individual explosions are not modeled. The 
effect of the supernova explosions is modeled on the sub-grid level 
as a blastwave propagating through the ISM. We track ISM metal 
enrichment from type II and type la supemovae as well as from 
AGB stars. 

The complete set of parameters used for the simulations dis- 
cussed in this Paper is listed in Table[T] Our fiducial run, has been 
studied extensively and we used it previously for results presented 
in R08ab. The baryonic particles in the fiducial run use a soften- 
ing length h a = 50 pc. We tested the softening dependence of our 
results with runs SI, S3 and S4 (h s — 25, 100 and 500 pc respec- 
tively). We also tested the effect of varying particle numbers in runs 
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Name 


N aas 


-^dark 


h s [kpc] 


fiducial 13 


10° 


10 6 


0.05 


SI 


10 6 


10 6 


0.025 


S3 


10 6 


10 6 


0.1 


S4 


10 6 


10 6 


0.5 


SDM 


10 6 


10 7 


0.05 


T2 C 


10 (i 


10 6 


0.05 


T3 C 


10 6 


10 6 


0.05 


Rl 


5 x 10 5 


5 x 10 5 


0.05 


R3 


2 x 10 6 


2 x 10 6 


0.05 


R4 


4 x 10 6 


4 x 10 6 


0.05 


Rl-T2 d 


5 x 10 5 


5 x 10 5 


0.05 



a All runs have M vir = 10 12 Mq, A = 
0.039 and c = 8.0. 

b This is our fiducial run that was presented 
previously in R08ab. 

c Runs T2 and T3 are identical to the fiducial 
run S2/R2 but use different random seeds when 
generating the initial conditions. 

d Run R1-T2 is identical to run Rl but initial- 
ized with a different random seed. 

Rl, R3 and R4 with 0.5 x 10 6 , 2 x 10 6 and 4 x 10 6 particles in 
each component respectively. We ran a further test of the effects of 
the mass resolution of the DM halo by running a simulation with 
10 times the number of DM particles (run SDM). All runs were 
carried out for 10 Gyr. 

The opening angle used for simplifying gravity calculations 
was 0.7 for all runs. The code uses a multistepping scheme with the 
condition that a particle's timestep At grav = ^(e/a) 1 / 2 , where 
r\ = 0.175, e is the gravitational softening length, and a is the 
acceleration. For the gas particles, the time step must also satisfy 

Atg as = r) coura nth/[(l + a) C + /3 fl m ax] , where Tlcourant = 0.4, 

h is the SPH smoothing length, a = 1 is the shear coefficient 
and = 2 is the viscosity coefficient. p, m ax is described in 
Wad slevetalj f2004). The SPH quantities were calculated on a ker- 
nel using 32 nearest neighbors. The cooling in all runs is calculated 
without taking into account the metal content of the gas. Additional 
simulations where we vary the physical parameters of the initial ha- 
los (mass, angular momentum) will be presented in Paper II of this 
series. 

The strength of these idealized models lies in the fact that 
there is no initial stellar component, so none of the properties of 
the disk are chosen a priori. As soon as the simulation begins, the 
gas is allowed to cool in the potential well of the DM, and as it 
reaches densities and temperatures conducive to star formation (as 
set by our star formation parameters) it spawns stars. Due to the 
spin, the gas naturally settles into a rotationally-supported disk in 
the center of the halo, giving rise to a galactic disk comprised of 
a mix of stars and gas. Thus, although all of the simulations pre- 
sented here are idealized in the sense that the disks evolve in iso- 
lation the simulations also depart sharply from most other work 
using idealized disks because we do not construct equilibrium disk 
models, a s is usually done in s tudies focusing on dynamical ef- 
fects (e.g. iDebattista et alj|2006l) . Instead, our models grow disks 
spontaneously without any biases regarding disk structure or stel- 
lar population properties. The spontaneous growth is particularly 
important for the development of self-limiting structure and subse- 
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quent secular evolution (e.g. ISellwood & Carlberdll984h . We are 
therefore able to follow the temporal evolution of stellar popula- 
tions and make direct comparisons between simulated properties 
and observed galaxies with unprecedented detail. 

Real disks of course do not form in such a simple way, but 
our simulations are meant only to mimic the formation of disks af- 
ter the last major merger. For massive galaxies (M ^ 10 M©) 
disk formation is not dominated by fila mentary accretion after z~2 
dKeres e t al. 20051: iBrooks et al.ll2009). Minor mergers may i nflu- 
ence the morphology of the disk (e.g. Kazantzidis et al. 200§), but 
the majority of the disk fuel comes from quiescent cooling of gas 
from the halo. 

By limiting the processes included in our models, we can per- 
form a cleaner analysis^ of the disk dynamics. Because of lower 
computational cost compared to fully-cosmological simulations, 
we are also able to attain higher resolution. At the fiducial resolu- 
tion, for example, the disks have ~ 2.5 million star particles by 10 
Gyr - typical hydrodynamic cosmological simulations of 1O 12 M0 
galaxies have a few times fewer disk particlefl However, by al- 
lowing our disks to form spontaneously without a pre-defined stel- 
lar component, we can follow the evolution of stellar populations 
self-consistently throughout the simulation. 



3 ANATOMY OF SPIRALS 

We begin by dissecting the spiral structure in our fiducial simula- 
tion to show that spirals are responsible for the migration of stars. 
For several decades, idealized TV-body simulations have produced 
self-limiting, transient spirals arising from various sources of in- 
stability, such as Poisson noise, d iscontinuities in the distribution 
function l Sellwood & Kahnll99lh . or stimulation by the central bar 
jDebattistaet^ll2006h . Studies of the evolution of transient spiral 
structure have largely been confined to the realm of carefully-tuned 
TV-body experiments, sometimes with an additi on of gas hydro- 
dynam ics or an ad-hoc cooling mechanism jSellwood & Carlberd 
11984) . Cooling is essential for the persistence of transient spirals 
because spirals are naturally self-destructive: when they form, they 
also very efficiently heat the disk. Controlled and carefully tuned 
TV-body experiments are vital for understanding the orbit response 
of the underlying disk to perturbing wave propagation, but do not 
adequately address how spontaneously recurring patterns may af- 
fect stars in a real disk. SB02 (their Figs. 10-12) provide a hint and 
show that an unconstrained collisionless simulation produces an ar- 
ray of spirals, identified by local maxima in the power spectrum of 
m=2 perturbations. In this section, we seek to expand on such anal- 
ysis by also considering the temporal evolution of patterns. 

In our simulations, the stellar disks are cooled naturally by on- 
going star formation. Stars are initially on kinematically cool orbits 
inherited from their parent gas particles, continuously infusing the 
disk with a population susceptible to supporting instabilities. The 
spiral structure manifests, as we show below, in a complicated ar- 
ray of patterns and pattern speeds that evolve as the disk grows. 
Figure Q] shows the Toomre Q = <jrk/3.36Y 1 G parameter as a 
function of radius at several different times, where or is the radial 



We analyze the simulations using our own IDL routines as well as the 
publicly available Python-based pynbody package: 
http :/ /code . google . com/p /pynbody / 

2 Cosmological simulations may have larger numbers of particles in the 
system, but most comprise the other galactic components like the bulge and 
halo and relatively few are found in the disk component. 
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R [kpc] 



Figure 1. Top: Toomre Q parameter as a function of radius at four different 
times in the fiducial ran. Bottom: Stellar surface density (black) and radial 
velocity dispersion (red). Line types correspond to the same times as in the 
upper panel. 



velocity dispersion of the stars, re is the epicyclic frequency, E is 
the surface density of stars, and G is the gravitational constant. The 
bottom panel shows £ and an separately. The recurring spirals are 
possible because the disk establishes a marginally stable state, with 
Q ~ 2 in the main part of the disk at all times. 

In order to study the spiral structure, we center the system on 
the potential minium and divide the particles into concentric equal- 
width radial bins. We then expand the stellar particle distribution in 
a Fourier series given by 



[ — imcf^Tn (r)] 



(2) 



where r is the radius, m is the pattern multiplicity, and <j> m {r) is 
the phase of the m-th mode at radius r. The complex coefficients 
c m (r) are given by 

1 N 

Cm(r) = MW^" 1 ^""^' (3) 



3=1 



where M(r) is the total mass in the radial bin, rrij is the particle 
mass, 4>j is the angle between the particle's position vector and the 
x-axis, and TV is the total number of particles in the bin satisfying 
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m=1 




± 



2 4 6 8 10 
time [Gyr] 

Figure 3. Global Fourier amplitudes of the stellar particles for m=l through 
m=5 Fourier components. 



r < rj < r + Sr. Strictly speaking, this type of decomposition 
will identify any type of azimuthal m-fold symmetry, but our as- 
sumption is that strong components identified in this way with low 
m will be due either to a bar or spirals. The bar can be separated 
from spirals by virtue of its constant phase as a function of radius; 
however, beyond the central few kpc we can reasonably expect this 
method to identify spirals. 

Figure [2] shows stellar density maps with contours of m=l 
through m=4 Fourier components, obtained using equations [2] and 
[3] demonstrating that our method can reliably identify and account 
for the underlying disk structure. It is also apparent from Figure [2] 
that there are often distinct spirals present at different radii and 
that each spiral often has a narrowly-defined peak overdensity. The 
overdensities we measure vary with time but the normalized ampli- 
tude of the m=2 mode A2 = \ct\ ^ 0.4 at all times for a ll runs, 
compatible with the observations of iRix & Z aritskv ( 1995). Over- 
densities at the level of a few per cent are difficult to see by eye 
in the stellar surface density map, but are easily identified with a 
Fourier decomposition. 

In a spontaneously-evolving disk, we expect patterns of all 
multiplicities to exist, and they all could perturb the orbits of stars. 
In Figure [3] we show the global Fourier amplitude as a function of 
time for m=l through m=5 Fourier components. The m=2 mode 
dominates at all times in this simulation. We therefore focus our 
dynamical analysis on the m=2 mode, though we cannot exclude 
that the other components also contribute some small fraction to 
the angular momentum redistribution in the disk. 

The greatest impact of spiral structure on the orbits of stars oc- 
curs at resonances, which are narrow regions in phase space satisfy- 
ing m(Q,p — Q p ) = Ik, where m is the pattern multiplicity, and 
k are the azimuthal and radial frequencies of the orbit and Sl p is the 
spiral's pattern speed. If m — 2, setting I to ±1 gives the inner and 
outer Lindblad resonances (ILR and OLR), and I = corresponds 
to the corotation resonance (CR). To identify resonances, we must 
therefore obtain fi p , which can be achieved using our Fourier de- 
composition. 



We can estimate the instantaneous Q p at radius r simply by 
performing numerical differentiation Q p = J^</!> m (r) or, equiva- 
lently, obtaining a linear fit to 4>m{r) as a function of time. This 
method is reliable if there is a single steady perturbation at the ra- 
dius in question, and studies of bars usually employ this method 
to recover the bar pat tern speed (e.g. lDebattista & SellwooJl200d : 
iDubinski et ail [2009). We find, however, that for our purposes 
this is insufficient because spirals of varying strengths and pattern 
speeds are present at all radii, and this method only ever identifies 
the most prominent pattern. 

To reliably identify multiple pattern speeds in the disk, we re- 
quire a different method. Ideally, the method would not only allow 
us to identify patterns at a given time interval, but it would also 
give us information about the evolution of these patterns with time. 
However, because the coefficients c m (r) are calculated at each out- 
put, they comprise a time series and we can obtain a further discrete 
Fourier transform of this series, given by 

s-i 

C*,m(r) = ^c j (r,m)w j e 2 * ijk/s k = 0,...,S-l, 
3=0 

< |Pressetal.ll 19921) . which yields the Fourier coefficients at discrete 
frequencies fife, and Cj (r, m) is the coefficient at timet = to+jAt, 
radius r and multiplicity m. S is the number of samples, in this case 
the number of outputs at which we evaluate the structure. We use a 
Gaussian window function 



to avoid high-frequency spectral leakage (using other window func- 
tions such as the Hanning window has no appreciable effect). The 
frequency sampling is determined by the length of the baseline and 
At, the time between samples (outputs), 

^ =27 W n fc = ' 1 -»2' 

where Qn v = fife=s/2 is the Nyquist frequency. We are here faced 
with a choice between spectral resolution and the ability to iden- 
tify "instantaneous" pattern speeds. If the baseline used is too long, 
the power spectrum will show a large number of patterns not all of 
which are important for the disk at any particular time. We find that 
at our timestep resolution of 10 Myr, using a 1 Gyr (100 timesteps) 
baseline gives satisfactory spectra, though the resolution is still 
rather coarse at ~ 3 km/s/kpc. This sampling rate gives a Nyquist 
frequency of £In v = 153 km/s/kpc for m — 2 structure. 

We can now construct a power spectrum at each radius given 

by 

P(Sl h ,r) = ± [\C k (r)\ 2 + \C S - k (r)\ 2 ] k = 1,2, ...f - 1, 

where W = S YLj-o w i ' s me normalization factor taking into ac- 
count the windowing function jPress et al.lll992h . We assume that 
the periodic changes will be due to the p hase, rather than amplitude 
variat ions, thus giving us pattern speeds dSellwood & Athanassoulal 
1986). The resulting combined power spectrum across the whole 
disk gives us information about the strengths of patterns as a func- 
tion of radius, and is shown in the Figure [4] The contours in the 
panels on the left show the frequency power spectrum as a function 
of radius, while the right panel shows the integrated power over the 
whole disk at each frequency. 

The pattern speeds in the given segment of the disk evolution 
can now be obtained by extracting the peaks from the radially inte- 
grated frequency spectrum. We measure the peaks by fitting Gaus- 
sians to the local maxima in the integrated spectrum. To reliably 
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Figure 2. Stellar density maps for the fiducial run at several times. Contours of overdensities reconstructed from the Fourier coefficients for m=l to m=4 are 
overlaid in red. Contours are drawn at the -50, -20, -5, 5, 20, and 50 per cent (under-)over-densities; negative contours are shown with dashed lines. We also 
require that each radial bin has at least 1000 particles to avoid noise that may arise in regions with low particle numbers. All panels are the same physical scale. 
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Figure 4. Left: Power spectrum of m=2 frequencies as a function of ra- 
dius in the disk, for the time interval 6-7 Gyr. The solid black line marks 
Q c , while the dashed black lines mark f2 e zfc/t/2. The three main spiral pat- 
terns are identified by colored horizontal lines. Right: The global frequency 
spectrum. Contour levels vary in logarithm from 5 X 10 — 10 to 1 X 10~ 6 . 



identify overlapping features in the spectrum we use an iterative 
procedure, which progressively removes the most prominent peaks 
identifying lesser peaks on subsequent passes. In this way, we iden- 
tify the three most important patterns. We use the Gaussian fit pa- 
rameters to estimate peak amplitudes and pattern speeds. The hor- 



izontal lines in the left panel of Figure [4] correspond to the pattern 
speeds obtained in this way and illustrate that this method can re- 
cover the major pattern speeds well. The solid and dashed black 
lines show Q c , the circular frequency, and fl c ± k/2 - ILR, CR, 
and OLR resonances exist for stars on nearly circular orbits where 
the pattern speed lines intersect one of these loci, from left to right 
respectively. 

Figure [4] reveals the richness of spiral structure in a sponta- 
neously evolving disk. The range of radii that may be influenced 
by resonances from any of these patterns spans essentially the en- 
tire disk. However, Figure|4]only tells us about the patterns present 
in a given time interval - we would like to know how these patterns 
evolve with time. 

To extract time information from our data, we repeat the above 
procedure at equally-spaced time intervals. We show such a result 
in the Figure [5] The dominant frequencies are shown as a func- 
tion of time, identified using the procedure outlined above, with 
individual point sizes reflecting the normalized amplitude of each 
perturbation. Each point on the plot uses a baseline of 1 Gyr of data 
to generate the power spectrum as discussed above, so that at each 
time t the patterns shown reflect disk evolution from t — 1 Gyr to 
t. Hence, adjacent points represent overlapping data and the point 
sizes reflect a time-averaged amplitude. The instantaneous ampli- 
tudes vary on much shorter timescales as the spirals flicker in and 
out of existence. 
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Figure 5. Perturbation frequencies in the disk as a function of time obtained 
via the windowed FFT (see text). Point sizes reflect the normalized ampli- 
tude of each perturbation. Color denotes the ratio of corotation radius to the 
break radius. 



It is striking that the perturbation at fl p ~ 20 — 
30 km/s/kpc dominates through most of the disk's lifetime. This 
perturbation is very steady except for a brief decrease at 5.5 Gyr 
and a subsequent strong growth peaking at 7.5 Gyr. The brief period 
of decreased power can also be identified in the panels of Figure|2] 
which show that between 5-6 Gyr the patterns become significantly 
disordered. Figure [5] also shows that perturbations with larger fre- 
quencies are significantly more short-lived. We can identify three 
separate strong perturbations with £l p > 50 km/s/kpc, the most 
persistent of which survives for ~ 2 Gyr. 

The time evolution of pattern speeds and amplitudes shown 
in Figure [5] provides a new look at the temporal evolution of the 
disk. We double-checked that the signal extracted using the Win- 
dowed FFT (WFFT) method is indeed real by performing an in- 
dependent analysis using a continuous wavelet transform (CWT), 
because it is specifically tailored to recover the time-evolution of 
frequencies in a signal. The CWT yields qualitatively identical re- 
sults to the WFFT - however, we found it less cumbersome to quan- 
tify the results (i.e. determine pattern speeds and amplitudes) using 
the WFFT. 

While Figure[5]shows the evolution of frequencies, it does not 
allow one to infer the amplitude variations of individual patterns 
because the WFFT is performed on 0.5 Gyr intervals. Therefore, 
while we do see variation in the sizes of the points, correspond- 
ing to amplitude changes, there is a possibility that we are hiding 
shorter time-scale fluctuations. 

We attempt to reconstruct the amplitude variation of individual 
patterns directly by using a band-pass filter in frequency space and 
then using the inverse Fourier transform to obtain the time series 
of Fourier coefficients representing the density perturbation. The 
band-pass is centered roughly on the desired frequency and is 10 
km/s/kpc wide. The choice of width is somewhat arbitrary, but we 
find that the patterns are reasonably well separated with this choice. 

Figure [6] shows the resulting amplitude variations for the three 
principal patterns found in Figure|4] The oscillatory amplitude tran- 
sience of spirals on Gyr timescales in our run is clearly visible. 
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Figure 6. Amplitude variation as a function of time for the three main pat- 
terns identified in Figure [4] The amplitudes are obtained by using a band- 
pass of 10 km/s/kpc centered on each frequency and at a radius correspond- 
ing to the corotation radius at 7 Gyr (i.e. identified using the power spectrum 
of Figure[4] Because the amplitudes of the faster patterns are quite different, 
we use the left y-axis for the 20-30 km/s/kpc pattern (red line), and the right 
y-axis for the other two (green and blue lines). 



Note that apart from the slowest pattern, the pattern speeds are also 
evolving on ~ Gyr timescales (see Figure [5). However, while the 
evolving pattern speeds mean that resonances associated with these 
patterns are sweeping through the disk, our system does not pro- 
duce a multitude of strong features of a wide range of pattern speeds 
at any given time. Thus, the transience in the system described here 
might be somewhat different from that found in other works. 



4 ANGULAR MOMENTUM EXCHANGES AT 
COROTATION IN THE SIMULATIONS 

In the previous section, we demonstrated that the disk in our simu- 
lation harbors a variety of spiral patterns of various pattern speeds. 
We now discuss how these patterns influence the orbits of stars. In 
particular, we try to understand whether the radial migration in our 
growing, unconstrained disk can be explained by the CR scattering 
described in SB02 alone, or whether we must invoke more complex 
mechanisms to explain the observed behavior. 

Figure UJ illustrates the diversity of individual stellar orbital 
histories in the disk. These orbits were chosen from a random set 
of particles that are found beyond 3 kpc at the end of the simula- 
tion, in order to illustrate the range of possible orbital histories. It 
is clear from these examples that stars may migrate rapidly while 
retaining a nearly circular orbit, and that a circular orbit does not 
imply a radially static history. In fact the instantaneous appearance 
of an orbit says very little about the star's history, as it is evidently 
possible to even circularize a significantly eccentric orbit (see, for 
example, the upper-right panel). 

In Figure [8] we show the changes in angular momentum as 
a function of initial particle angular momentum^ For the bottom 
panel, the perturbations dominating the disk during this time are 
the same ones we identified in Figure[4] Note that the baseline used 
to obtain the frequency power spectrum is 1 Gyr, but the timescale 
on which we consider Aj z vs. j z is 0.5 Gyr, where j z — [r x v] 



3 Throughout the paper we use the notation that specific angular momen- 
tum is lower-case j and total angular momentum is capital J 
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Figure 7. Some examples of orbits from the fiducial simulation. Orbits are colored according to time to make the temporal evolution more apparent, with blue 
corresponding to early and red to late times. Stars can migrate inwards and outwards very rapidly without substantially increasing their eccentricities. 



and Ajz — jz {ti ) — j z (ti). If we use a 1 Gyr baseline (which cor- 
responds to > 4 rotation periods at 5 Gyr and 8 kpc) we find that the 
signature of spiral scattering in the Aj z vs. j z plane becomes too 
smeared to be useful. On the other hand, a 0.5 Gyr baseline is in- 
sufficient to reliably identify pattern speeds in the power spectrum, 
so we are forced to accept this small inconsistency. The higher fre- 
quency perturbations may evolve considerably during ~ 1 Gyr, but 
the resonance loci are not a particularly steep function of 7? (Fig- 
ure so we can still obtain approximate locations of resonances. 
In Figure [6] we can see that during 6.5-7 Gyr we may expect the 
disk to be largely dominated by the growth of the 45 km/s/kpc pat- 
tern. However, during the 5.3-5.5 Gyr time interval, we see that the 
two faster patterns are growing at the same time. It is also appar- 
ent from Figure[5]that around this time the faster pattern dominates 
the overall perturbation amplitude in the main part of the disk. We 
therefore choose this second time interval to look for nonlinear ef- 
fects steirHTung_froirij^sM multiple perturbations 
(e.g. lOuilledl2003klMinchev & Famae\ll2010h . 

In both cases, the dominant pattern causes the characteris- 
tic anti-symmetric shape in Aj z vs j z associated with CR orbit 
swapping (SB02). The few dominant patterns we identify cannot 
explain the full structure seen in either panel, however. In Fig- 
ure [4] we can see that some amount of power is present in sub- 
tler patterns not picked up by our automated algorithm, such as at 
~ 30 km/s/kpc and ~ 60 km/s/kpc. Some of the structure in Fig- 
ure[4]must therefore also be due to those patterns, which may be in 
the process of growing or fading at either end of the time interval 



we consider. Another possibility is that patterns of higher multi- 
plicity contribute to some of the angular momentum exchange. The 
prominent feature seen in the very inner part of the disk in ther bot- 
tom panel is a short-lived weak bar instability. Note the scale of 
Aj z : in both panels the changes caused by the dominant pattern 
approach 50% of j z . The changes are especially drastic in the top 
panel, given the short time interval we are considering. 

4.1 Corotation Crossing and Chaos 

A steady spiral perturbation does not result in angular momentum 
exchange at the CR because stars are trapped on horseshoe orbits. 
If the spiral is transient, however, the trapping does not occur and 
stars can traverse from one side of the CR to another. We look for 
such horseshoe orbits during the time interval 6.5-7 Gyr among 
the particles with |Aj z j > 200 (roughly the top 10 percent of mi- 
grators) around the CR with the 43 km/s/kpc pattern shown in the 
bottom panel of Figure [8] We plot the orbits of a random subset of 
these particles in a frame corotating with the spiral in columns 1 
and 3 of Figure [9] The orbits are colored according to the relative 
strength of the 43 km/s/kpc spiral at any given time, shown in the 
bottom right panel. The radius as a function of time is shown as 
a black line in columns 2 and 4. In the corotating frame, a parti- 
cle crossing the CR also reverses direction. All particles cross the 
CR at almost the same time, near the peak of the spiral amplitude, 
but the orbits are otherwise largely unperturbed. Furthermore, all 
inward and outward migrations are respectively correlated in az- 
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Figure 8. Distributions of Aj z given initial j z for all particles in the sim- 
ulation between 5.3-5.5 Gyr (top panel) and 6.5-7.0 Gyr (bottom panel). 
The corotation resonances of the main perturbations found at each time are 
marked by vertical lines, with the colors in the top panel corresponding to 
the colors in Figure [4] The logarithmically-spaced contour levels are the 
same for both panels and correspond to mass density. 



imuth. For example, see the outward migrators in panels 1, 2 and 
3 and the inward migrators in panels 6 and 7. This hints at the CR 
being responsible for the migration because one expects the parti- 
cles trailing the spiral arm to preferentially be pulled outward and 
vice versa for those particles leading the arm. 

In Figure [10] we demonstrate that this is indeed the case by 
plotting the density of inward and outward migrators (in red and 
blue respectively - darker colors indicate higher density) for the 
5.3-5.5 Gyr time interval (top panel) and 6.5-6.7 Gyr time inter- 
val (bottom panel). Spiral overdensity, reconstructed from m = 2 
through m — 4 Fourier components, is shown in black contours. 
The first interval is chosen because it is near the time when two 
of the patterns are peaking simultaneously, while during the later 
time interval only the 43 km/s/kpc spiral peaks. Note that we did 
not select these particles specifically to isolate ones interacting with 
either spiral. The only selection criterion is that in the given time 
interval they are in the top 5 percent of migrators. However, essen- 
tially all of the migration in both time intervals is occuring around 
the CRs with most of the inward and outward migrators respec- 
tively leading and trailing the dominant spiral. Low densities of 
migrators extend to other regions of the disk where other spirals 
play a role, but the majority of migration is taking place around the 



43 km/s/kpc spiral, as expected from Figure [8] The dashed green 
circle corresponds to the CR of the 43 km/s/kpc spiral and passes 
directly through the centers of the highest density regions of migra- 
tors. 

The top panel, showing the migrators for the time interval 
where two of the faster patterns are peaking simultaneously shows 
very similar behavior to the bottom panel where only one pattern is 
dominating. If the overlap of resonances was causing chaotic evolu- 
tion, it should manifest itself during such a time interval (the domi- 
nant patterns are ~ 65 and 30 km/s/kpc). Around a radius of 5 kpc 
in the top panel, the density of both inward and outward migrators 
is saturated, implying that the migration is even more dominated 
by a single pattern in this time interval. The physical locations of 
the migrators are strongly clustered around the individual overden- 
sity peaks. In this interval, there is also a weak oval feature in the 
center, driving some angular momentum exchange in the inner part 
the disk (see the top panel of Figure [8). The dashed green circle in 
this panel corresponds to the CR of the 65 km/s/kpc spiral, which 
is clearly where most of the angular momentum exchange is oc- 
curing. Thus it is evident that even when multiple strong patterns 
are present near their amplitude peaks, the majority of angular mo- 
mentum exchange occurs at very predictable locations, i.e. on ei- 
ther side of overdensity peaks and at their CR. The migration in 
this earlier time interval is stronger, hinting at the possibility that 
different patterns may help feed each other's CR zones. However, 
as we show in the following paragraphs, we are unable to find clear 
indications of non-linear or chaotic evolution. 

In columns 2 and 4 of Figure [9] we show the continuous 
wavelet transform (CWT) power scalogram of the x component 
of the orbits shown in columns 1 and 3. We construct the power 
scalogram by using the Moiiet-Grossmann wavelet whose basis 
"mother" wavelet is given by 



-t 2 /2o 



See lDaubechiesI dl990h for a thor ough t heoretical introduction to 



Torren ce & Compol ( 1 19981) and iNener et all dl999h 



wavelets and 

for a practical guide. We use the approach described in the lat- 
ter two for calculating our wavelet power scalograms. The y-axis 
in the spectrograms in Figure [9] represents wavelet scale, which 
is proportional to the inverse of frequency. We can therefore fol- 
low the time evolution of dominant frequencies by inspecting the 
most prominent ridges in the scalogram. Note that we plot the log- 
arithm of the power, so many of the smaller features are orders 
of magnitude weaker than the obvious dominant ridges. Chaotic 
systems can be identified with this method because the ridges in 
the scal ogram show si gnificantly irregular curvatu re and disconti- 
nuities (Chandr e et alj2003l ; lGemmeke et alj2008h . 

Chand re et al. 1 20031) show that in chaotic systems, ridges in 
the scalogram become shorter, bent, and disordered. In our case, we 
find in general mostly horizontal features extending across much 
of the entire time interval during which migration takes place. As 
the stars cross corotation and migrate in radius, the dominant fre- 
quency has to change and this is reflected in the bend apparent in 
the main ridge in all the panels. However, there is no evidence of 
significantly scattered and bent ridges (see for example Fig ure 4 in 
Gemmeke et al. 2008 and Figure 13 of lChandre etal]|2003h . 

The lack of obvious evidence for chaos is particularly impor- 
tant for particles that migrate significantly in a short amount of 
time. To investigate this further, we looked at the orbits of the two 
particles in the top row of Figure [7] which migrate ~ 5 kpc in 
0.5 Gyr. These two particles are chosen for their extreme rates of 
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Figure 10. Density of particles migrating significantly in the time intervals 
5.3-5.5 Gye (top) and 6.5-6.7 Gyr (bottom). The sole selection criterion is 
that the particles are in the top 5 percent in terms of their | Ajr' z | in the given 
time interval. The blue and the red colors show the outward and inward 
migrators respectively. The black contours show the surface density of stars 
reconstructed from the m = 2 through m = 4 Fourier components. The 
dashed green circles mark the CR of the dominant pattern in the given time 
interval (65 km/s/kpc in the top panel and 45 km/s/kpc in the bottom). The 
direction of rotation is counter-clockwise. 

migration, where we may also expect to detect hints of nonlinear- 
effects. In the left columns of the top two rows of Figure [TJJ we 
show the the orbits corotating with the faster, ~ 65 km/s/kpc pat- 
tern, and in the middle panels the orbits are in the frame of the 
~ 40 km/s/kpc spiral. As in Figure(9] the orbits are colored accord- 
ing to the relative strength of the two spirals, whose amplitudes are 
shown in the bottom row. The radius as a function of time for each 
particle is shown as a black line in the rightmost column. Both par- 
ticles cross the CR at almost the same time near the peak of the 
spiral amplitude, similar to the particles in Figure [9] The two par- 
ticles are first taken to the middle of the disk by the faster pattern 
and happen to be deposited near the corotation resonance and near 
the peak of the spiral pattern dominating there. They are therefore 
quickly taken further outward by this second, slower spiral. Thus, 
their rapid migration is simply due to two successive corotation res- 
onances. As before, the wavelet transform power spectrum shown 
in the rightmost column of Figure QT] shows little evidence for sig- 
nificantly chaotic evolution. 

We look further for signs of nonlinearity by calculating the 
change in the Jacobi integral for a subset of the migrators. The Ja- 



cobi integral, Ej = E — Q p j z is conserved in the rotatin g frame 
of a single, steady perturbation ( Binnev & Tremainell2008l) . In our 
system we do not expect Ej to be exactly conserved anywhere since 
multiple patterns are always present in the disk, and the disk poten- 
tial itself is constantly changing due to accretion of fresh gas and 
star formation. However, if single spirals are mostly responsible for 
instantaneous changes in j z of individual stars, we would expect 
that for short time intervals the distribution of AE/ Aj z — Q p to 
be peaked and centered at zero for particles in the vicinity of each 
of the major CR regions. We select particles for these distributions 
by determining the j z ,GR, i.e. the angular momentum of circular 
orbits at each CR, and then selecting the top 10 percent of migra- 
tors with initial j z within 10 percent of j z ,cn- In essence, we are 
selecting narrow strips of particles on either side of each CR line 
shown in Figure [8] Note that we determined the pattern speeds in 
the 6.5-6.7 case separately from those determined in Figure [4] and 
there is a slight discrepancy due to pattern speed evolution over Gyr 
timescales. 

The distributions for both time intervals are shown in Fig- 
ure [12] All distributions are centered near zero, within our un- 
certainties of pattern speed measurements of a few km/s/kpc. The 
distributions for the slower perturbations dominating in the outer 
disk are very strongly peaked, as expected if most of the particles 
are only being affected by a single spiral. The distributions for the 
faster perturbations are broadened for several reasons. First, the spi- 
ral amplitudes are variable which broadens their effective resonant 
frequencies. Second, in the inner disk other weaker patterns may 
be affecting the orbits. Still, it is striking that even in the 5.3-5.5 
Gyr time interval, which was chosen specifically because multiple 
strong patterns are peaking at the same time (see also bottom panel 
of Figure [Tot, the distributions do not appear qualitatively any dif- 
ferent than for the quieter 6.5-6.7 Gyr time interval. The higher 
density of migrators seen in Figure [10] could be due to several fac- 
tors. First, the disk density is significantly higher at the CR of the 
60 km/s/kpc spiral. The angular momentum exchange could fur- 
ther be facilitated by the presence of other perturbations, essentially 
feeding the CR of the dominant spiral. Nevertheless, even if the ef- 
ficiency of migration is boosted due to the presence of other pat- 
terns, the fact that the distributions shown in Figure[l2]are centered 
and peaked at zero implies that the bulk of the angular momentum 
exchange still proceeds due to the CR of the individual spirals and 
that chaotic orbital evolution caused by the overlap of several reso- 
nances may not be the dominant cause of migration in this system. 
Note that the shift to the left from zero for the red distribution in 
the bottom panel of Figure[l2]is just a few km/s/kpc and therefore 
well within the measurement errors for f2 p . 

Similarly, we examine the evolution of Ej for the two par- 
ticles shown in Figure [TT] that undergo very large migration over 
a short time interval. In the bottom right panel we show the frac- 
tional changes in the Jacobi integral for these two particles. The 
fractional change for the faster pattern is measured with respect to 
the beginning of the time interval, whereas for the slower pattern 
we measure it with respect to the end of the interval (this is because 
initially the particles interact with the faster pattern and only later 
with the slower pattern). Between 5-5.6 Gyr, Ej in the frame of 
the 65 km/s/kpc pattern varies by < 1%. In the same time interval, 
the energy E (shown in red) and angular momentum (not shown 
for brevity) of these two particles change by 20% and 50% respec- 
tively. Between 5.6-6 Gyr, Ej in the frame of the ~ 40 km/s/kpc 
pattern (blue lines) is again roughly constant, while energy and an- 
gular momentum continue to evolve and change by 30% and a fac- 
tor of 2.5 by the end of the time interval. 
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Figure 9. Columns 1 and 3: A selection of randomly selected orbits for a subsample of top 10 percent of migrators on either side of the CR for the 43 km/s/kpc 
pattern in the time interval 6.5-7 Gyr. The orbits are in the frame corotating with the spiral. The color corresponds to the relative spiral amplitude, as shown 
in the bottom right panel. All boxes measure 20 kpc in x and y. The direction of rotation is counter-clockwise. Note that the time of CR crossing is easily 
identified as the orbit reverses direction. All CR crossings happen around the time of the spiral amplitude peak and are spatially correlated. Columns 2 and 
4: The wavelet power spectrum of the rr-component of the orbits is shown in color contours. Note that most of the features are horizontal with few bends and 
vary smoothly, indicating lack of significant chaotic evolution. The solid black lines show the radius R in the xy plane as a function of time. See the bottom 
left panels for axis scales. 



4.2 Orbital Circularity and Migration 



We now focus on the issue of orbital circularity and heating in CR 
migration. In Figure [13] we show that the amount of angular mo- 
mentum exchange is directly related to the circularity of an orbit. 
This is a crucial feature of corotation scattering, i.e. particles on the 
most circular orbits are also the most susceptible to having their or- 
bits drastically altered. It follows that the largest changes in angular 
momentum are those that are essentially kinematically untraceable, 
because corotation scatt ering does not increase the ran dom energy 
of an orbit appreciably jLvnden-Bell & Kalnaisll 19721 SB02). Par- 
ticles are selected in the same way as for Figure [12] (for brevity, 
we only show the middle pattern that causes the largest amount 
of mixing). From the rotation curve and the midplane potential, 
we calculate the theoretical circular orbit locus in the (j, E) plane. 
Based on this locus we calculate the maximum allowable angular 
momentum, j c {E), of each particle based on its energy and express 
the circularity of its orbit as jz/jc(E). 

We calculate the mass weighted distribution of |Aj 2 | as a 



function of x = jz/jc(E) which can be expressed as 

N 

f( X ) = J2\ A jz,i\ m i S ( X - X i), ( 4 ) 
i=l 

where the subscript i indicates individual particle quantities, m is 
the mass, and TV is the total number of particles. The corresponding 
cumulative distribution function is given by 

= (5) 

and shown in Figure [T3] Note that even for particles on fairly ec- 
centric orbits, changes in j z are possible, but those are most likely 
occurring at the Lindblad resonances. However, the contribution of 
those stars to the overall changes in j z is miniscule. The particles 
with jz/jc{E) > 0.95 account for over 50% of the angular mo- 
mentum exchange. This follows since | Aj z | is largest for those par- 
ticles, but also simply due to the fact that the disk is kinematically 
cool and particles on mostly circular orbits are also most abundant. 
However, the mass distribution (shown by the dashed line in the 
same panel) is much less peaked than the angular momentum dis- 
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Figure 11. Columns 1 and 2: Similar to Figure[9]but for the two orbits shown in the top row of Figure |7] The orbits are in the frame corotating with the 
inner spiral (column 1) and the mid-disk spiral (column 2) in the time interval 5-6 Gyr, when they happen to migrate substantially. The color corresponds to 
the relative spiral amplitude, as shown in the bottom row. All boxes measure 20 kpc in x and y. The direction of rotation is counter-clockwise. Note that the 
time of CR crossing is easily identified as the orbit reverses direction. Column 3: The top two rows show the wavelet power spectrum of the ^-component of 
the orbits in color contours. The solid black lines show the radius R in the xy plane as a function of time. The bottom panel in column 3 shows the fractional 
time evolution of the Jacobi constant, Ej in the frame of the two spirals (green and blue lines correspond to the 65 and 39 km/s/kpc patterns respectively) and 
energy E (red lines). Due to different ranges of values, the left- and right-hand y-axes are used for Ej and E respectively. Particles 1 and 2 are shown in solid 
and dashed lines. When the particles are in the corotation region of a given pattern, Ej is conserved to ~ 1%, while E evolves by 30% (though not shown, 
angular momentum changes by a factor of 2.5) 



tribution, confirming that particles on circular orbits are most im- 
portant for the exchange of angular momentum. 

Apart from being most efficient at relocating particles on the 
most circular orbits, SB02 also showed that CR scattering also does 
not appreciably heat the disk. In Figure [14] we show a probability 
density distribution of A[j z / j c (E)] given AR over the same time 
interval studied in Figures[8]and[T2] Positive values of A[j z / j c (E)] 
indicate an increase in circularity or a "cooling" of the orbit, while 
negative values correspond to heating. The majority of the parti- 
cles that migrate outwards suffer very minor heating, while the 
inward-migrating particles heat slightly more. The particles mov- 
ing the farthest outward also get heated the least on average - this 
is the key feature of CR scattering because it allows particles to ex- 
perience multiple scatterings, thereby allowing for potentially very 
large changes in radius during their lifetimes. Interestingly, roughly 
10% of the particles that migrate outward have their orbits cooled 
by the spiral. 

Note from the top panel of Figure [8] that the most obvious 
and dramatic feature in Aj z vs. j z occurs due to a resonance with 



the 43 km/s/kpc pattern that does not appear dominant in Figure[5] 
The ~20 km/s/kpc spiral that dominates the power spectrum on the 
other hand causes relatively little mixing. Determining the dom- 
inance of patterns from Figure [4] is somewhat misleading because 
the power spectrum is constructed using the normalized Fourier co- 
efficients. Thus, although the slowest pattern appears to dominate, 
it peaks in the far outer disk where the density is low and therefore 
the perturbation does not involve much mass. The middle pattern 
instead is also at an optimal point in the disk. The dashed line in 
Figure[T]shows the Toomre Q parameter at 8 Gyr, which is near the 
interval we have been analyzing. In the inner disk where the fast 
pattern CR occurs, Toomre Q approaches values of 4 and above - 
similarly in the outer disk, where the slow pattern CR occurs Q in- 
creases rapidly and the disk density is very low (the break occurs at 
~ 8 kpc but the CR is at ~ 11 kpc). For the middle pattern CR, Q 
is relatively low while the disk density on the other hand is still rea- 
sonably high. Therefore the CR is well-populated by kinematically 
cool stars and the middle pattern can achieve the most mixing. 
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Figure 12. Distributions of AE/AJ Z — fl p for particles with the largest 
changes in J z near the corotation resonances of patterns identified from the 
Fourier analysis. The top and bottom panels show the 5.3-5.5 and 6.5-6.7 
Gyr time intervals respectively. 
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Figure 13. Cumulative distribution of jz/jc{E) weighted by Aj z \ - see 
equations [4] and \5\ for a description of F[x), The dashed line shows the 
cumulative distribution of j z /j c (E) . 
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Figure 14. Probability distribution of the change in the ratio of jz/jc{E) 
for the same particles used to construct the green distribution in the top 
panel of Figure ll2l and in Figure [T3l interacting with the 42 km/s/kpc pattern 
over the 6.5 - 7 Gyr time interval. The inner (outer) white contours enclose 
50% (75%) of the mass in each AR bin. Values are normalized such that in 
the vertical direction, ^2 PdA = 1, where dA is the area of the bin. 



5 NUMERICAL TESTS 

Spiral arms are amplified disturbances in disks, but the source of 
the seed perturbation is not well understood. It is often assumed 
that in is olated disks, the seed is noise in the density distribu- 
tion (e.g. Goldre ich & Lvnden-Belllll965l) . In simulations like the 
ones presented here, the natural source of noise is the finite par- 
ticle numbers, which are in general several orders of magnitude 
smaller than the number of stars in real galactic disks, though 
in real galaxies giant molecular clouds are a similar source of 
"noise". Numerical resolution studies abound in the literature, but 
most often attention has been given to the requirements of col- 
lisionless cos mological simulati ons to resolve dark matter sub- 
structure (e.g. lMoore et al.| [T998). In the cases where disk secular 
evolution has been addressed specifically, most resolution studies 
have addressed the resonant couplings between dark matter ha- 
los and bars, b ut again most often only collisionless simulations 
were used (e.g. Debattista & Sellwood 2000; Valenzuela & Klvpi 
20031; iHollev-Bockelmann et alj 120051 ; IWeinberg & Katzl bOO 
Dubinski et alJl2009T) . Spiral structure and in particular SPH sim 



IIS 

7: 



ulations in isolation ar e addressed less often, though recently 
Christensen et al. I J2010h presented a global disk resolution study 
focusing mostly on the resolution dependence of sub-grid star 
formation and feedback prescriptions. Unfortunately, their high- 
resolution runs use approximately the same resolution as the fidu- 
cial runs here, because their work was meant as a guide for cosmo- 
logical simulations where the state-of-the-art uses comparable or 
slightly lower resolution. Their results therefore cannot be used to 
ascertain the validity of our results. 

In this section, we discuss several numerical tests to explore 
the robustness and variability of the spiral activity presented in the 
preceding sections. The main concern here is that the resulting spi- 
rals may be dependent on the Poisson particle noise for their gener- 
ation and subsequent evolution. We performed runs with 0.5, 2, and 
4 times the fiducial particle number of 10 6 particles per component 
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in the initial conditions. Softened gravity could be particularly rele- 
vant to disks because it can potentially set the relevant perturbation 
scales. We therefore ran several simulations with different choices 
for the softening length. A possible source of numerical heating of 
the disk are massive dark matter particles, so we also ran a simu- 
lation with an order of magnitude more dark matter particles. We 
discuss each of these in turn below. 

5.1 Stochasticity 

ISellwood & Debattistal d2009h examined the effects of stochasticity 
on the growth and pattern speed evolution of bars. They found that 
small perturbations in the initial conditions could lead to divergent 
behavior, independent of the code used for the integration. Spirals 
result from amplified disturbances and are intrinsically sensitive to 
stochastic effects. It is therefore not possible to expect the pattern 
speed and amplitude evolution to match exactly among the differ- 
ent runs. To get some sense of the natural range of behavior due to 
stochasticity, we ran two simulations with identical numerical pa- 
rameters and initial conditions, but we altered the random seed used 
in the generation of initial particle positions sampled from the dis- 
tribution function. The general disk properties that result from 10 
Gyr of evolution are very similar, with inner scale lengths ranging 
from 3.1 to3.3kpcQ. 

In Figure [15] we show the evolution of the m=2 Fourier am- 
plitude at four different radii in the disk for all runs. We show 
smoothed m=2 time series because otherwise the rapid oscillations 
of the amplitudes make it difficult to discern their overall evolu- 
tion. The colors represent the amplitude at different radii. Bar/oval 
growth can be identified in this representation whenever the black 
and blue lines grow together. The middle and right panels of the top 
row show the stochasticity test runs (fiducial, T2 and T3) - appar- 
ently the timing of the growth of spirals is very different between 
the three runs with different random seeds. Run T2 does not grow a 
central oval at all (the black and blue lines showing 2 and 4 kpc are 
almost completely featureless), whereas the fiducial run and run T3 
both have episodes of bar/oval formation. In the fiducial run, this 
occurs at ~ 3, 4, and 6 Gyr, whereas for the T3 run it occurs at 3 
and 6 Gyr. In Figure [16] we show the pattern speed evolution de- 
rived in the same way as in Figure[5]with the panels corresponding 
to panels in Figure[T5] Very similar pattern speeds occur in all three 
disks, as shown by the top row of Figure[T6] For example, the same 
outer (slow) pattern can be seen in the three power spectra, starting 
at around 30 km/s/kpc and slowly decaying to 20 km/s/kpc by 10 
Gyr. The pattern speeds are qualitatively similar in all runs and pre- 
sumably depend upon the overall disk structure, which is relatively 
robust against stochastic effects. The amplitude and timing of the 
spirals, on the other hand, may vary considerably. 

5.2 Softening 

The second row of Figure Q3] shows the time evolution of the m=2 
perturbations at four different radii for the runs with different values 
of the softening parameter. This suggests that the simulations with 
h s ^ 100 pc yield qualitatively similar spiral structure, though SI 
(leftmost panel) harbors somewhat more damped spirals, especially 

4 These values were obtained using fits to surface density profiles. If we 
instead fit midplane volume-density profiles, the resulting scale lengths are 
~ 2.5 kpc, in agreement with ljuric et alj 120081) values obtained for the 
Milky Way from SDSS data. 



mid-disk (green line). Run S4 is clearly wildly different, develop- 
ing a large and persistent bar very early (in subsequent discussions 
we ignore this run). This comparison shows significant variance in 
the timing of spiral activity, but the amplitudes at different radii 
for runs SI and S3 appear very similar to the fiducial run (upper 
left panel). The second row of Figure [T6l shows the corresponding 
pattern speed evolution and confirms that for runs SI and S3, the 
pattern speeds supported by the disk are very similar to each other 
and to the fiducial run. 

5.3 Particle Number 

The third row of Figures Q3] and [Tfj] show the m=2 amplitude and 
pattern speed evolution for runs Rl, R3, and R4 (0.5, 2, and 4 times 
the fiducial particle number respectively). Run Rl stands out in this 
comparison, as the structure that develops is significantly weaker. 
We investigated whether this is a manifestation of stochasticity and 
ran a second simulation where the initial conditions were generated 
using a different random seed (run R1-T2). The results are shown in 
the right panel of the bottom row - the evolution of this experiment 
is essentially identical to Rl, suggesting that the weaker structure 
is not a stochastic effect. 

The star formation rates are the same in all four runs. As a 
result, the disks are comprised of proportionally similar amounts 
of stars - roughly 2, 1,4, and 8 million for the fiducial, Rl, R3, and 
R4 runs respectively. 

In runs R3 and R4, the structure is much more similar to the 
fiducial run. The timing of perturbation growth is slightly different, 
but the frequencies are essentially the same for all four runs. The 
consistency of these results affirms that the runs in our simulation 
suite have sufficient numbers of particles to adequately model the 
disk dynamics. Ideally, we would be able to perform simulations 
with still higher particle numbers to truly test for convergence, but 
due to computational cost such simulations were not feasible at this 
time. 

The ratio of CR to the break radius (i.e. the color of the points) 
for the ~ 30 km/s/kpc pattern in run R4 appears different from 
t ~ 6.5 Gyr onwards. These differences are not particularly drastic, 
however, and the relative locations of the CR are consistent. The 
speeds of the main patterns are similar to the runs in the top row 
(fiducial, Tl, T2). We also see an even slower pattern, probably 
due to increased particle numbers, which make its detection in the 
far outer disk possible. 

Finally, because dark matter particles are constantly bombard- 
ing the disk, we explored the possibility that their perturbations 
may influence the generation of spirals. In the bottom-left panel of 
Figure [T5l we show the m=2 amplitude for run SDM, which was 
initialized with 10 times as many DM particles. Other properties of 
the run were kept the same as in the fiducial run. Spiral amplitudes 
are very similar to the fiducial run, although the innermost pattern 
seems to be missing. In the corresponding panel in Figure [16] the 
frequencies as a function of time for the run SDM appear very sim- 
ilar to the fiducial run. At most times, the dominant frequencies 
present in both disks are essentially the same. 

While we have achieved convergence with increasing mass 
resolution, the discrepancy between our lowest resolution run, Rl, 
and the rest of the suite is puzzling. We tested for the effect of 
two-body relaxation by increasing the softening and by increasing 
the number of particles in the DM halo. We also tested for poten- 
tial integration issues stemming from different softenings between 
the baryons and the DM (note that in all of our runs the DM and 
baryons use different softening, and this is standard practice in such 
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Figure 15. Fourier amplitudes as a function of time for the entire suite. The fiducial run is shown in upper left, followed by the stochasticity tests in the same 
row. Second row shows runs with different softenings; the third row shows runs with different particle numbers. Bottom row shows runs SDM and R1-T2. 
Different color lines represent the amplitude at different radii - black, blue, green, and red correspond to 2, 4, 8, 12 kpc respectively. 



simulations). Making the baryon softening equal to the DM soften- 
ing was the only experiment that yielded results more in line with 
other runs. 



5.4 Overall Comparison 

In the preceding subsections we investigated the detailed evolu- 
tion of spirals in our simulation suite, but how do these simula- 
tions compare in their global properties? In Figure[T7]we show the 
mean m=2 amplitude, (A2), calculated over the duration of the en- 



tire simulation, and its standard deviation represented by the error 
bars. Stochasticity and choice of softening minimally impact the 
global disk evolution, though all of the S-series runs have slightly 
smaller (A2), especially S4. A true outlier is Rl, with lower (A2) 
than the rest of the suite. 

The properties of spirals are important in our understanding 
of radial migration because of their effect on the distribution of 
stars in the disk. We have seen in Figures [15] and [16] that while 
the range of admissible pattern speeds is not sensitive to choices of 
numerical parameters, the timing of the spirals can differ substan- 
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Figure 16. Pattern speed evolution as a function of time. Colors and point sizes are as in Figure[5] Panels correspond to panels from Figure [T31 



tially. Depending on the configuration of the disk, this could affect 
the migration rates and ultimately, the predictions we extract from 
these models for studies of disk stellar populations. It is difficult to 
quantitatively assess the spiral structure, though we have attempted 
to do so in Figure [T7] and found no appreciable differences in the 
suite apart from run Rl. However, we can quantitatively analyze 
the properties of resulting stellar populations in a given region of 
the disk. 

A natural region to examine is the model analog to the solar 
neighborhood. In the top panels of Figure[T8]we show the distribu- 
tions of formation radii for stars that are found within 7 < i?[kpc] < 
9 at the end of the simulation. This choice is particularly important 



becaus e we have addressed t his same region in our previous works 
(R08b, Loebma n et al.ll2oTTh . The bottom panels show the corre- 
sponding cumulative distributions. The solid black line in all panels 
corresponds to the fiducial run - the left, middle, and right columns 
show the stochasticity, softening, and resolution tests respectively. 

The leftmost panels show that stochasticity has little overall 
effect on the cumulative distributions of Rjorm - at 50% the dif- 
ference is ^ 0.5 kpc. Nevertheless, the distribution in the top panel 
shows that the fiducial run may even over-produce the in-situ stars. 

In the remaining cumulative distributions, the overall variance 
at 50% does not exceed that of the stochasticity tests. Two notable 
cases are apparent - the run SI and run SDM. Run SI, is more 
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Figure 17. The mean m=2 amplitude over the whole simulation for all of 
the runs in our suite. The error bars indicate the standard deviation. Points 
are offset horizontally for readability. 



heavily dominated by in-situ stars. This agrees with the fact that 
this run develops weaker spiral structure in the final Gyr. However, 
the appearance of the larger peak of in-situ stars for the run SI 
may also simply be another manifestation of stochasticity - when 
we recreate the same distributions 2 Gyr earlier, the run S 1 follows 
almost exactly the fiducial distribution. 

For run SDM, the peak of the Rform distribution is actually 
shifted away from the solar neighborhood. This implies even more 
drastic mixing, and is indeed also apparent when we repeat the ex- 
periment at 8 Gyr (as mentioned in the previous paragraph). Re- 
gardless of these subtle variations, the fact that > 50% of solar 
neighborhood stars originated in other parts of the disk remains ro- 
bust. 

Contrary to what we might expect for run Rl given Figures[T5l 
[l6]and[T7] where we found that it has weaker spiral structure on av- 
erage from other runs, we find here only a small difference in the 
solar neighborhood population. We would expect the resultant mi- 
gration to be much less than in the other runs. We therefore studied 
the spiral structure in this run in more detail and found that while 
the disk does not support strong spirals, it is permeated with weaker 
and highly transient features. We find that the transience is much 
more pronounced than in the other runs. Using a 1 Gyr baseline 
power spectrum, as in Figure [4] we find essentially a continuous 
spectrum of patterns emerges, compared with few distinct features. 
This is apparent also in Figure [76] where many weak patterns at 
a variety of pattern speeds are identified. The pattern speeds shift 
on short timescales, allowing even the weak structure to relatively 
efficiently redistribute the stars. 



6 CONCLUSIONS 

Radial migration is rapidly becoming recognized as a critical pro- 
cess in disk galaxies because of its wide-ranging implications. 
If this process has been important for the Milky Way, it most 
likely played a key role in determining the mix of stars in the 
solar n eig hborhood today {Sellwood & Binnevll2002l; Roskar et al.l 



2008b al: ISchonrich & Binnevl l2009dlbl; [Loebman et alj 1201 it 
Minchev & Famaevl l2010t |Bird et alj |2012| ; iLee et al.l 1201 lh . On 
a broader scale, it has likely influenced the stellar population 



gradients measured in the Milky Way disk and o ther galaxies 
(iBoissier & Prantzoj 120001: iMacArthur et all |2004 IRoskar et aT 



2008 



Williams et alj 120091: iGogarten et alj l201Ct IVlaiic et al 



201ll ;[Munoz M ateos et al.ll201lh , and contributed substantially to 



the stellar den sity in the outermost disk reg i ons dBarker et a l. 200 



dejongetal. 2007; Roskar et al. 2008b; Azzollini et al. 2008 



iBakos et alj|200a IRoskar et al.ll 2010; Yoachim et alj|201fl 2012, 
in press; lBarker et alj|201 ltlAlberts et al 1201 lh . 

In this Paper, we have investigated in detail the origin of radial 
migration by analyzing the spontaneously-forming spiral structure 
and the resultant resonant angular momentum exchange. We found 
that the spiral structure is transient in amplitude, but appears to 
support only a few discrete pattern speeds at any given time. This 
means that some stars can be tossed from the CR of one pattern 
to another, resulting in large changes in radius on relatively short 
timescales. Still, it is important to remember that extreme migra- 
tions of many kpc over the course of a star's lifetime are not the 
norm, they comprise the tail of the distribution. This can be seen in 
Figure [T8l- although ~50% of the stars do come from elsewhere, 
this also means that the near-majority have formed in-situ. The sit- 
uation changes with increasing radius, because the in-situ star for- 
mation decreases - thus, the tail of the distribution from the in- 
ner regions makes up the majority of the stellar population at large 
radii. 

We demonstrated in Section|4]that the largest angular momen- 
tum exchanges occur at the corotation of important m=2 spirals. An 
important aspect of our result is the confirmation that the largest 
angular momentum exchange happens for particles on the most cir- 
cular orbits - and that these particles do not get heated by the spi- 
rals while they migrate (SB02). This is a crucial aspect of the CR 
migration mechanism because it means that the process is not self- 
limiting. Instead, it can continue especially for the particles migrat- 
ing the most, allowing for very large changes in radius for some of 
the stars, but without betraying their journey by anomalous kine- 
matics. We have also searched for signs of chaotic orbital evolution 
in the vicinity of resonances, but found indication that stellar orbits 
remain rather regular as they migrate radially. We have found this 
to be the case even for very large migrations, and have shown that 
such migrations are possible in a short amount of time if a particle 
passes directly from the CR of one spiral to another (see Figures [9] 
andE}. 

Our results show that the redistribution of stars in MW-type 
disks on very short timescales is inevitable if transient spiral struc- 
ture is present. Even in the case of run Rl, where the mid-disk spi- 
rals appear to be numerically suppressed, a large fraction of stars 
ending up in the solar neighborhood originated in the interior of 
the disk. Therefore, our findings suggest that the effects of recur- 
ring, spontaneous spiral structure are a key component of disk evo- 
lution that models simply must include if they wish to make pre- 
dictions about kinematic and chemical properties of stellar popula- 
tions. Cosmological simulations, which fail to form disks that can 
support spiral structure may be missing critical aspects of disk evo- 
lution and therefore the detailed properties of resulting disk stellar 
populations must be used with care. On th e other hand, the inclu- 
sion of substructure (e.g. iBird et al.lr2012T) is important not only 
because mergers can heat and disrupt the disk, but also because 
they may trigger transient structure. We hope that it should soon 
be possible to use state-of-the-art cos mological simulations (e.g. 
iGuedes et all 1201 it lAgertz et al]|201lh for detailed studies of ra- 
dial migration. 
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